Meine Gedanken zur Workshopgestaltung
- Möglichst viel Aktivität bei den Teilnehmer*innen
- Weniger Vorturnen durch mich, mehr selbst ausprobieren (und erstmal scheitern 😃)
- Ich zeige zunächst alles in , jeder benutzt danach seine Lieblingssoftware - ich helfe bei Fragen zu JASP/jamovi/SPSS/STATA/etc. so gut ich kann …
- Möglichst viel Interaktion
- Fragen zu Begriffen gerne direkt stellen
- Wünsche wiederholter/weiterer Erklärungen oder Elaborationen gerne direkt äußern
- Gerne Peer-to-Peer Interaktion in Break-Out-Rooms
- Gemeinsames Interpretieren der Lehrbuchtexte (in Sitzung 2)
- Möglichst differenziertes Arbeiten
- Gestufte Lösungshilfen
- Break-Out-Rooms für “aktiv Modellierende” und “passiv Nachvollziehende”
- Kognitiv aktivierende Inputs die auf verschiedenen Verständnisebenen rezipiert werden können
- Materials
- Datensätze aus .html downloadbar
- Für -Userinnen gibt es hier das komplette RStudio Projekt: https://github.com/sammerk/Universaltool_Regressionsanalyse_II_-Mehrebenenmodelle-.git
- Die Dokumentation kann im Browser unter https://bit.ly/merk035 abgerufen werden. Alle Veränderungen werden dort auch sichtbar
- Real World Data (führt oft zu Softwareproblemen etc. aber das entspricht der Realität)
Worked Out Example: Interpretation multipler Regressionsmodelle
Datensatz 1: Kid IQ
Der Datensatz Kid IQ stammt aus einem der empfohlenen Lehrbücher(Gelman & Hill, 2007). Ursprünglich stammt er aus der National Longitudinal Survey of Youth. Wir werden die folgenden Variablen Nutzen:
- kid_score = Rohwert des Kidnes in einem Intelligenztest
- mom_hs = Dummyvariable Highschoolabhschluss (1 = Highschoolabschluss, 0 = kein Highschoolabschluss)
- mom_iq = IQ der Mutter
- mom_age = Alter der Mutter
- mom_work =
- 1: mother did not work in first three years of child’s life
- 2: mother worked in second or third year of child’s life
- 3: mother worked part-time in first year of child’s life
- 4: mother worked full-time in first year of child’s life
Import der Daten
library(tidyverse)
data_kidiq <- read_delim("https://raw.githubusercontent.com/sammerk/did_data/master/kidiq.csv", delim = ";")
Diejenigen, die nicht nutzen wollen, können den Datensatz hier im SPSS-Format herunterladen.
Modelle mit einem metrischen Prädiktor
Parametrisierung
Zunächst soll der kid_score mit der metrischen Variable
mom_iq prädiziert werden. In
erfolgt das mit der Syntax
mod01 <- lm(kid_score ~ mom_iq, data = data_kidiq)
mod01
##
## Call:
## lm(formula = kid_score ~ mom_iq, data = data_kidiq)
##
## Coefficients:
## (Intercept) mom_iq
## 25.80 0.61
Grafisch kann dieses Modell wie folgt repräsentiert werden:
library(hrbrthemes)
data_kidiq %>%
ggplot(., aes(mom_iq, kid_score)) +
geom_point() +
stat_smooth(method = "lm", se = F) +
theme_ipsum()
## `geom_smooth()` using formula 'y ~ x'
Und die formale Beschreibung lautet:
\[ \operatorname{\widehat{kid\_score}} = 25.8 + 0.61(\operatorname{mom\_iq}) \]
Da der kid_score offensichtlich nicht einer
Standardskalierung (z-Werte, IQ-Werte, t-Werte, …) entspricht, kann die
Stärke des Effekts nur anhand einer Standardisierung bewertet
werden.
mod02 <- lm(scale(kid_score) ~ scale(mom_iq), data = data_kidiq)
mod02
##
## Call:
## lm(formula = scale(kid_score) ~ scale(mom_iq), data = data_kidiq)
##
## Coefficients:
## (Intercept) scale(mom_iq)
## -5.068e-16 4.483e-01
Dies leisten auch packages wie {sjPlot} oder
{stargazer} die darüber hinaus auch noch mehrere Modelle
vergleichend darstellen können.
library(sjPlot)
tab_model(mod01, show.std = T, show.ci = F, show.p = F)
| kid_score | ||
|---|---|---|
| Predictors | Estimates | std. Beta |
| (Intercept) | 25.80 | -0.00 |
| mom iq | 0.61 | 0.45 |
| Observations | 434 | |
| R2 / R2 adjusted | 0.201 / 0.199 | |
Nach Cohen (1988) gilt dabei:
- \(\beta \approx .1 \Rightarrow \text{kleiner Effekt}\)
- \(\beta \approx .3 \Rightarrow \text{moderater Effekt}\)
- \(\beta \approx .5 \Rightarrow \text{starker Effekt}\)
Inferenzstatistik
Unser \(\beta = .45\) beschreibt
lediglich die Steigung einer Geraden, wenn man diese nach dem Kriterium
der kleinsten Quadrate auf dem Datensatz kid_iq optimiert.
Möchte man Schlussfolgerungen über den datengenerierenden Mechanismus
anstellen benötigt man Inferenzstatistik.
p-Werte kann man etwa anhand der summary()-Funktion oder
der tab_model() Funktion erhalten.
summary(mod01)
##
## Call:
## lm(formula = kid_score ~ mom_iq, data = data_kidiq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -56.753 -12.074 2.217 11.710 47.691
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.79978 5.91741 4.36 1.63e-05 ***
## mom_iq 0.60997 0.05852 10.42 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.27 on 432 degrees of freedom
## Multiple R-squared: 0.201, Adjusted R-squared: 0.1991
## F-statistic: 108.6 on 1 and 432 DF, p-value: < 2.2e-16
Dabei sind zwei Inferenzstatistiken enthalten:
- Jeder Koeffizient \(\beta_i\) wird gegen die 0 getestet.
- Das Gesamtmodell wird gegen \(R^2 = 0\) getestet.
JZS-Bayes-Faktoren sind via {BayesFactor}
ermittelbar:
library(BayesFactor)
lmBF(formula = kid_score ~ mom_iq, data = data_kidiq)
## Bayes factor analysis
## --------------
## [1] mom_iq : 4.516119e+19 ±0%
##
## Against denominator:
## Intercept only
## ---
## Bayes factor type: BFlinearModel, JZS
Diagnostik der Voraussetzungen für die Inferenzstatistik
Diese Inferenzstatistiken sind nur unter Annahmen über den datengenerierenden Mechanismus berechenbar. Sind diese Annahmen verletzt sind die Inferenzstatistiken (mehr oder weniger stark) verzerrt. Daher nimmt die Diagnostik der Annahmen eine zentrale Bedeutung ein. Angenommen werden muss (mindestens)
- Linearität der Beziehung
- Normalverteilung der Residuen
- Homoskedastizität der Residuen
- Unabhängigkeit der Residuen
Eine gute Grundlage für die Diagnostik bietet das
{performance}-package:
library(performance)
check_model(mod01)
Interpretation
Regressionsmodelle beschreiben immer nur bedingte Erwartungswerte von Datenpunkten (“für eine Gruppe von Müttern deren IQ im Schnitt X ist, ist ein kid_score von X am wahrscheinlichsten”). Woher diese Wahrscheinlichkeit rührt und wie gut mit dieser eine kausale Relationierung von X und Y gerechtfertigt werden kann hängt maßgeblich vom Design der Studie ab!
Modelle mit einem dichotomen Prädiktor (Dummyvariable)
Parametrisierung
Die Regression ist insofern ein recht universelles Modellierungstool,
als dass es auch die Aufnahme dichotomer Prädiktoren erlaubt. Bspw. kann
man mit der Variable mom_hs der kid_score
prädizieren:
ggplot(data_kidiq, aes(mom_hs, kid_score)) +
geom_point(alpha = .3) +
stat_smooth(method = "lm", se = F) +
theme_ipsum()
## `geom_smooth()` using formula 'y ~ x'
Die resultierenden Koeffizienten stellen die arithmetischen Mittelwerte der Gruppen dar und der p-Wert von \(\beta_1\) tested die \(H_0: EW(Gruppe_1) = EW(Gruppe_2)\).
mod03 <- lm(kid_score ~ mom_hs, data = data_kidiq)
mod03
##
## Call:
## lm(formula = kid_score ~ mom_hs, data = data_kidiq)
##
## Coefficients:
## (Intercept) mom_hs
## 77.55 11.77
Wobei \[ \operatorname{kid\_score} = \beta_{0} + \beta_{1}(\operatorname{mom\_hs}) + \epsilon \]
Diagnostik der Voraussetzungen für die Inferenzstatistik
Da Regressionsmodelle mit Dummyvariablen eben auch Regressionsmodelle sind 😃 gilt es auch bei diesen dieselben Voraussetzungen zu prüfen:
check_model(mod03)
Modelle mit mehreren Prädiktoren (multiple Regression)
Parametrisierung
In Regressionsmodelle können problemlos mehrere Prädiktoren aufgenommen werden.
mod04 <- lm(kid_score ~ mom_iq + mom_hs, data = data_kidiq)
summary(mod04)
##
## Call:
## lm(formula = kid_score ~ mom_iq + mom_hs, data = data_kidiq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -52.873 -12.663 2.404 11.356 49.545
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.73154 5.87521 4.380 1.49e-05 ***
## mom_iq 0.56391 0.06057 9.309 < 2e-16 ***
## mom_hs 5.95012 2.21181 2.690 0.00742 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.14 on 431 degrees of freedom
## Multiple R-squared: 0.2141, Adjusted R-squared: 0.2105
## F-statistic: 58.72 on 2 and 431 DF, p-value: < 2.2e-16
Grafisch kann dies dann im entsprechenden n-dimensionalen Raum repräsentiert werden (siehe Folien aus Universaltool Regressionsanalyse I). In der formalen Repräsentation wird schlicht das Polynom um einen Summanden erweitert: \[ \operatorname{kid\_score} = \beta_{0} + \beta_{1}(\operatorname{mom\_iq}) + \beta_{2}(\operatorname{mom\_hs}) + \epsilon \]
Interpretation
Besonders interessant an multiplen Regressionsmodellen ist, dass bestimmte kausale Mechanismen falsifiziert werden können. Dazu werden dann meist die Ergebnisse der einfachen Modelle (mit jeweils einem Prädiktor) mit den Parameterschätzungen im multiplen Modell verglichen:
tab_model(mod01, mod03, mod04, show.ci = F)
| kid_score | kid_score | kid_score | ||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | p | Estimates | p | Estimates | p |
| (Intercept) | 25.80 | <0.001 | 77.55 | <0.001 | 25.73 | <0.001 |
| mom iq | 0.61 | <0.001 | 0.56 | <0.001 | ||
| mom hs | 11.77 | <0.001 | 5.95 | 0.007 | ||
| Observations | 434 | 434 | 434 | |||
| R2 / R2 adjusted | 0.201 / 0.199 | 0.056 / 0.054 | 0.214 / 0.210 | |||
So können Annahmen über Supressor- oder Konfounderkonstellationen falsifiziert werden. Hier etwa wird klar, dass die Daten mit dem Mechanismus A nicht im Einklang stehen, aber mit dem Mechanismus B und C.
Diagnostik der Voraussetzungen für die Inferenzstatistik
Auch Modelle der multiplen Regression unterliegen den gleichen 4
Annahmen. Hinzu kommt (je nach Methode der Schätzung), dass keine
Kolinearität vorliegen sollte. check_model() nimmt dies
Annahme direkt in die Prüfung mit auf, wenn das Argument der Funktion
ein Modell mit mehreren Prädiktoren darstellt.
check_model(mod04)
Typischerweise werden dann wieder p-Werte für die Hypothesen \(\beta_i = 0\) und $ R^2 = 0$
berechnet.
Bayes Faktoren beziehen sich jedoch immer auf das inkrementelle \(R^2\), also die Frage: Klärt ein Prädiktor
oder eine Gruppe von Prädiktoren zuätzliche Varianz auf?
# Test gegen intercept only modell
lmBF(formula = kid_score ~ mom_iq + mom_hs, data = data_kidiq)
## Warning: data coerced from tibble to data frame
## Bayes factor analysis
## --------------
## [1] mom_iq + mom_hs : 1.508721e+20 ±0%
##
## Against denominator:
## Intercept only
## ---
## Bayes factor type: BFlinearModel, JZS
# Test gegen Modell mit mom_hs als Prädiktor
lmBF(formula = kid_score ~ mom_iq + mom_hs, data = data_kidiq)/
lmBF(formula = kid_score ~ mom_hs, data = data_kidiq)
## Warning: data coerced from tibble to data frame
## Warning: data coerced from tibble to data frame
## Bayes factor analysis
## --------------
## [1] mom_iq + mom_hs : 7.578673e+15 ±0%
##
## Against denominator:
## kid_score ~ mom_hs
## ---
## Bayes factor type: BFlinearModel, JZS
# Test gegen Modell mit mom_iq als Prädiktor
lmBF(formula = kid_score ~ mom_iq + mom_hs, data = data_kidiq)/
lmBF(formula = kid_score ~ mom_iq, data = data_kidiq)
## Warning: data coerced from tibble to data frame
## Warning: data coerced from tibble to data frame
## Bayes factor analysis
## --------------
## [1] mom_iq + mom_hs : 3.340747 ±0%
##
## Against denominator:
## kid_score ~ mom_iq
## ---
## Bayes factor type: BFlinearModel, JZS
Modelle mit Interaktionseffekt eines metrischen und eines dichotomen Prädiktors
Parametrisierung
Interessiert man sich dafür, wie sich die Effekte je nach Ausprägung einer weiteren nominalen Variablen unterscheiden, kann man Modelle mit Interaktionen eines metrischen und eines dichotomen Prädiktors spezifizieren.
ggplot(data_kidiq, aes(mom_iq, kid_score,
color = as.factor(mom_hs),
group = as.factor(mom_hs))) +
geom_point() +
stat_smooth(method = "lm") +
theme_ipsum()
## `geom_smooth()` using formula 'y ~ x'
Die entsprechenden Koeffizienten des Modells mit Interaktionseffekt sind
jedoch wesentlich schwieriger zu interpretieren:
mod05 <- lm(kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq, data = data_kidiq)
mod05
##
## Call:
## lm(formula = kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq, data = data_kidiq)
##
## Coefficients:
## (Intercept) mom_hs mom_iq mom_hs:mom_iq
## -11.4820 51.2682 0.9689 -0.4843
\[ \operatorname{kid\_score} = \alpha + \beta_{1}(\operatorname{mom\_hs}) + \beta_{2}(\operatorname{mom\_iq}) + \beta_{3}(\operatorname{mom\_hs} \times \operatorname{mom\_iq}) + \epsilon \]
Leichter wird dies, wenn man den kontinuierlichen Prädiktor zentriert oder standardisiert.
library(equatiomatic)
data_kidiq <- data_kidiq %>%
mutate(mom_iq_centered = as.numeric(scale(mom_iq, center = T, scale = F)),
mom_iq_zstand = as.numeric(scale(mom_iq, center = T, scale = T)),
mom_age_zstand = as.numeric(scale(mom_age, center = T, scale = T)))
mod06 <- lm(kid_score ~ mom_hs + mom_iq_zstand + mom_hs:mom_iq_zstand, data = data_kidiq)
mod06
##
## Call:
## lm(formula = kid_score ~ mom_hs + mom_iq_zstand + mom_hs:mom_iq_zstand,
## data = data_kidiq)
##
## Coefficients:
## (Intercept) mom_hs mom_iq_zstand
## 85.407 2.841 14.533
## mom_hs:mom_iq_zstand
## -7.264
ggplot(data_kidiq, aes(mom_iq_zstand, kid_score,
color = as.factor(mom_hs),
group = as.factor(mom_hs))) +
geom_point() +
stat_smooth(method = "lm") +
theme_ipsum()
## `geom_smooth()` using formula 'y ~ x'
extract_eq(mod06, use_coefs = T, terms_per_line = 2, wrap = T)
\[ \begin{aligned} \operatorname{\widehat{kid\_score}} &= 85.41 + 2.84(\operatorname{mom\_hs})\ + \\ &\quad 14.53(\operatorname{mom\_iq\_zstand}) - 7.26(\operatorname{mom\_hs} \times \operatorname{mom\_iq\_zstand}) \end{aligned} \]
Diagnostik der Voraussetzungen für die Inferenzstatistik
Auch hier gelten dieselben Annahmen, die wieder gut mit
check_model() evaluierbar sind:
check_model(mod06)
Modelle mit Interaktionseffekt zweier metrischer Prädiktoren
Betrachtet man Modelle mit Interaktionen zwischen zwischen zwei metrischen Variablen, sind diese am besten interpretierbar, wenn beide zentriert sind.
mod07 <- lm(kid_score ~ mom_iq_zstand*mom_age_zstand, data = data_kidiq)
mod07
##
## Call:
## lm(formula = kid_score ~ mom_iq_zstand * mom_age_zstand, data = data_kidiq)
##
## Coefficients:
## (Intercept) mom_iq_zstand
## 86.8761 9.2097
## mom_age_zstand mom_iq_zstand:mom_age_zstand
## 1.1404 -0.8633
extract_eq(mod07, use_coefs = T, wrap = T)
\[ \begin{aligned} \operatorname{\widehat{kid\_score}} &= 86.88 + 9.21(\operatorname{mom\_iq\_zstand}) + 1.14(\operatorname{mom\_age\_zstand}) - 0.86(\operatorname{mom\_iq\_zstand} \times \operatorname{mom\_age\_zstand}) \end{aligned} \] In diesen Modellen können dann die sogenannten Haupteffekte, also die \(\beta_i\) der einfachen Summanden ohne multiplikative Terme (hier \(\beta_1\) und \(\beta_2\)) als die Effekte der Variable \(i\) gesehen werden, wenn die andere Variable eine durchschnittliche Ausprägung (nach Zentrierung = 0) hat.
Richtig hübsch geplottet wird das per default mit der Funktion
plot_model()
plot_model(mod07, type = "int") +
theme_ipsum()
Übung: Multiple Regression mit Dummyvariablen
Datensatz 2: Effekte der Klassengröße und Wertzuschreibung
Das Student Teacher Achievement Ratio (STAR) Projekt ist eines der
größten Experimente zu Effekten der Klassengrößenreduktion auf die
Schüler*innenleistung. Die Klassen wurden randomisiert drei Bedingungen
zugewiesen: Große Klasse, kleine Klasse, große Klasse mit
Hilfslehrkraft. Für die dritte Klasse findet sich diese Information in
der Variable g3classtype (1 = SMALL CLASS, 2 = REGULAR
CLASS, 3 = REGULAR + AIDE CLASS). Außerdem wurde über einen
Lehrerfragebogen die Wertzuschreibung der Schülerinnen und Schüler
erfasst. Im Datensatz ist die Variable mit dem Kürzel
g4ptvalu zu finden. Ein Beispielitem lautet This
student thinks that school is important mit den
Antwortmöglichkeiten 1 = Never, 2, 3 = Sometimes, 4, 5 = Always. Die
Mathematikleistung ist für die dritte und vierte Klasse in den Variablen
g3tmathss und g4tmathss enkodiert.
Die Daten können wie folgt heruntergeladen werden (ein .sav-fiel gibt es